
load('C:/Users/michael/Dropbox/NFLIS/statelevel.Rdata')
load("C:/Users/michael/Dropbox/NFLIS/itemlevel.Rdata")

library(data.table)
library(ggplot2)
library(geofacet)

################################ start itemlevel

itemlevel = data.table(itemlevel)
itemlevel$YEAR = as.numeric(substr(itemlevel$submitdate, 7, 10))

itemlevel$anyfentanyl = as.numeric(grepl("Fentanyl", paste(itemlevel$drug1new,
              itemlevel$drug2new, itemlevel$drug3new, itemlevel$drug4new,
              itemlevel$drug5new, itemlevel$drug6new, itemlevel$drug7new,
             itemlevel$drug8new))
)

itemlevel$anyamph = as.numeric(grepl("Amphetamine", paste(itemlevel$drug1new,
                                                           itemlevel$drug2new, itemlevel$drug3new, itemlevel$drug4new,
                                                           itemlevel$drug5new, itemlevel$drug6new, itemlevel$drug7new,
                                                           itemlevel$drug8new))
)

itemlevel$anycocaine = as.numeric(grepl("Cocaine", paste(itemlevel$drug1new,
                                                          itemlevel$drug2new, itemlevel$drug3new, itemlevel$drug4new,
                                                          itemlevel$drug5new, itemlevel$drug6new, itemlevel$drug7new,
                                                          itemlevel$drug8new))
)

itemlevel = itemlevel[(anyfentanyl + anyamph + anycocaine) > 0]
itemlevel$cocaine_fentanyl = as.numeric((itemlevel$anyfentanyl + itemlevel$anycocaine) ==
                                          2)
itemlevel$amph_fentanyl = as.numeric((itemlevel$anyfentanyl + itemlevel$anyamph) ==
                                          2)



item_stateyear = itemlevel[, list(cocaine_fentanyl=sum(cocaine_fentanyl),
                                  amph_fentanyl=sum(amph_fentanyl),
                                  cocaine_total=sum(anycocaine),
                                  amph_total=sum(anyamph)
                                  ),
                           by=list(state, YEAR)]

############## END ITEMLEVEL

statelevel = data.table(statelevel)
statelevel$YEAR = statelevel$CaseReceived_Year

statelevel$anyfentanyl = as.numeric(grepl("Fentanyl", paste(statelevel$drug1new,
                        statelevel$drug2new, statelevel$drug3new, statelevel$drug4new))
)

statelevel$anyamph = as.numeric(grepl("Amphetamine", paste(statelevel$drug1new,
                                  statelevel$drug2new, statelevel$drug3new, statelevel$drug4new
))
)

statelevel$anycocaine = as.numeric(grepl("Cocaine", paste(statelevel$drug1new,
                                  statelevel$drug2new, statelevel$drug3new, statelevel$drug4new
))
)

statelevel = statelevel[(anyfentanyl + anyamph + anycocaine) > 0]
statelevel$cocaine_fentanyl = as.numeric((statelevel$anyfentanyl + statelevel$anycocaine) ==
                                          2)
statelevel$amph_fentanyl = as.numeric((statelevel$anyfentanyl + statelevel$anyamph) ==
                                       2)


statelevel$state = statelevel$State
state_stateyear = statelevel[, list(cocaine_fentanyl=sum(cocaine_fentanyl),
                                  amph_fentanyl=sum(amph_fentanyl),
                                  cocaine_total=sum(anycocaine),
                                  amph_total=sum(anyamph)
), by=list(state, YEAR)]

#####################

## combine
combo_stateyear = rbind(item_stateyear, state_stateyear)

# re-iterate sum
combo_stateyear =  combo_stateyear[, list(cocaine_fentanyl=sum(cocaine_fentanyl),
                                     amph_fentanyl=sum(amph_fentanyl),
                                     cocaine_total=sum(cocaine_total),
                                     amph_total=sum(amph_total)
), by=list(state, YEAR)]


### plot
plot_data = melt(combo_stateyear,
                 measure.vars = c("amph_total", "cocaine_fentanyl",
                "amph_fentanyl", "cocaine_total", "amph_total"))
plot_data = plot_data[variable %in% c("cocaine_fentanyl", "amph_fentanyl"),]
ggplot(plot_data, aes(YEAR, value, color = variable, linetype=variable)) +
  geom_point() + ggtitle("Fentanyl-Stimulant Seizures (2011-2016)") +
  geom_line() + facet_geo(~state) + xlab("") + 
  scale_x_continuous(breaks = seq(2012, 2016, 2)) +
  ylab("Annual seizures") + theme_bw() + labs(color="Fentanyl with",
                                              linetype = "Fentanyl with") +
  scale_color_discrete(label = c("Cocaine", "Methamphetaine")) +
  scale_linetype_discrete(label = c("Cocaine", "Methamphetaine")) +
  theme(plot.title = element_text(size=28), axis.text = element_text(size=8))
ggsave("fent_amph_cocaine_overlay.png", dpi=300, width=12, height=8, units="in")

ggplot(plot_data[variable=="amph_fentanyl",], aes(YEAR, value, color = variable)) +
  geom_point() + ggtitle("Fentanyl-Methamphetamine Seizures (2011-2016)") +
  geom_line() + facet_geo(~state) + xlab("") + 
  ylab("Annual seizures") + scale_x_continuous(breaks = seq(2012, 2016, 2)) +
  theme_bw(base_size = 12) + theme(legend.position = "",
  plot.title = element_text(size=28), axis.text = element_text(size=8))
ggsave("fent_amph.png", dpi=300, width=12, height=8, units="in")

ggplot(plot_data[variable=="cocaine_fentanyl",], aes(YEAR, value, color = variable)) +
  geom_point() + ggtitle("Fentanyl-Cocaine Seizures (2011-2016)") +
  geom_line() + facet_geo(~state) + xlab("") + 
  ylab("Annual seizures") + scale_x_continuous(breaks = seq(2012, 2016, 2)) +
  theme_bw(base_size = 12) + theme(legend.position = "",
  plot.title = element_text(size=28), axis.text = element_text(size=8))
ggsave("fent_cocaine.png", dpi=300, width=12, height=8, units="in")

# stimulant
combo_stateyear$stim = combo_stateyear$cocaine_fentanyl + combo_stateyear$amph_fentanyl
ggplot(combo_stateyear, aes(YEAR, stim)) +
  geom_point() + ggtitle("Fentanyl-Stimulant Seizures (2011-2016)") +
  geom_line() + facet_geo(~state) + xlab("") + 
  ylab("Annual seizures") + scale_x_continuous(breaks = seq(2012, 2016, 2)) +
  theme_bw(base_size = 18) + theme(legend.position = "",
  plot.title = element_text(size=28), axis.text = element_text(size=8))
ggsave("fent_combinedstim.png", dpi=300, width=12, height=8, units="in")

